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Abstract 

We investigate second order linear wave equations in periodic media, 
aiming at the derivation of effective equations in M. n , n > 1. Standard 
homogenization theory provides, for the limit of a small periodicity length 
e > 0, an effective second order wave equation that describes solutions 
on time intervals [0, T]. In order to approximate solutions on large time 
intervals [0,Te~ 2 ], one has to use a dispersive, higher order wave equation. 
In this work, we provide a well-posed, weakly dispersive effective equation, 
and an estimate for errors between the solution of the original heterogeneous 
problem and the solution of the dispersive wave equation. We use Bloch- 
wave analysis to identify a family of relevant limit models and introduce an 
approach to select a well-posed effective model. The analytical results are 
confirmed and illustrated by numerical tests. 

Keywords: homogenization, wave equation, weakly dispersive model, 
Bloch-wave expansion 
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1 Introduction 

The wave equation describes wave propagation in very different applications, rang- 
ing from elastic waves to electro-magnetic waves. In the two mentioned applica- 
tions, it is of interest to describe waves in periodic media where the period is much 
smaller than the lengthscale of the wave. The most fundamental questions regard 
the effective wave speed and a possible dispersive behavior due to heterogeneities. 

We concentrate on the simplest model, the second order wave equation in 
divergence form. For notational convenience we omit here a positive coefficient in 
the acceleration term and study, for x G M. n , the wave equation 

}2„ e 



<9 t V(x, t) = V • (a e (x)Vu e (x, t)) . (1.1) 



Technische Universitat Dortmund, Fakultat fur Mathematik, Vogclpothsweg 87, D-44227 
Dortmund, Germany. 



1 



2 



T. Dohnal, A. Lamacz, B. Schweizer 



The medium is characterized by the coefficient function a £ : R n — > (0, oo). Since 
we are interested in periodic media with a small periodicity length-scale e > 0, 
we assume that a £ (x) = ciy{x/e), where ay : MJ 1 — > R has the periodicity of the 
cube Y := (— vr,7r) n C M n . The wave equation is complemented with the initial 
condition 

u £ (x,0) = f(x), d t u £ (x,0) = 0. (1.2) 
We consider initial data / G L 2 (R n ) n L x (K n ) of the form 

f(x) = j^^F (k)e^dk, (1.3) 

and assume that the (bounded) function F : M n — > C has a compact support 
if CC R". In this sense, / has a finite support in Fourier space. 

The fundamental question of homogenization theory is: For small e > 0, 
can the solution u £ be approximated by a solution of an equation with constant 
coefficients? The answer is affirmative: There exists an effective coefficient matrix 
A G M nx?1 , computable from ay, such that the following holds: on an arbitrary 
time interval [0, T], if we define w : R n x [0, T] — > R as the solution of 

d 2 t w{x,t) = V • (AVw(x,t)) , w(x,0) = f(x), d t w{x,0) = 0, (1.4) 

there holds u £ — > w as e — > 0. For the result and function spaces see e.g. [4]. 

We are interested in a refinement of this result. Our aim is to investigate the 
behavior of solutions u £ of (1.1) for large times. To formulate a precise result, 
we start from a positive number To > 0, and investigate solutions for all t G 
[0,To£~ 2 ]. It is well-known that the homogenized equation (1.4) cannot provide 
an approximation of u £ on the interval [0,To£ -2 ]. Instead, we need a dispersive 
equation to approximate u £ . 

Main result. In addition to the coefficient matrix A, we will define matri- 
ces E and F. All matrices are computable from the coefficient ay(.) with the 
help of an eigenvalue problem on the periodicity cell Y. The constant coefficient 
matrices define linear spatial differential operators, the two second order opera- 
tors AD 2 = . Aijdidj and ED 2 = z~2i j Eijdidj, and the fourth order operator 
FD A = Y^ijmi Fijmididjdmdi . We formulate a weakly dispersive equation 

d 2 w £ = AD 2 w £ + e 2 ED 2 d 2 w £ - e 2 FD 4 w £ . (1.5) 

As initial conditions we use once more w £ (x, 0) = f(x) and d t w £ (x, 0) = 0. Equa- 
tion (1.5) is of fourth order in the spatial variables, but it contains additionally an 
operator which involves second spatial and second time derivatives. The operator 
contains the small parameter e > explicitly; formally, for e = 0, we recover the 
homogenized equation (1.4). 

Our main result shows that the weakly dispersive equation (1.5) provides, at 
the lowest order, an approximation of the original equation for large times. 
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Theorem 1.1. Let e = Si — > be a sequence of positive numbers, n > 1 the di- 
mension. Let the heterogeneous medium be given by a Y -periodic positive function 
ay : M n -}Ro/ class C 2 , symmetric under reflections Xj < — > —Xj, and symmet- 
ric under coordinate exchanges Xj < — > x^, see (2.23). Let initial data be given by 
f G L l (M. n ), such that the Fourier transform F has compact support, see (1.3). 
We assume that the Bloch- expansion (2.8) of the initial data f is convergent in 
H l (W l ). 

We use the coefficient matrices A and C defined in (2.19), and E and F as 
defined in Lemma 3.1. Then the following holds: 

1. Well-posedness Equation (1.5) with initial condition (1.2) has a unique 
solution w £ for all positive times (see Theorem 3.1 for function spaces). 

2. Error estimate Let w e be the solution of (1.5), and let u e be the solution 
of (1.1) for the same initial condition (1.2). Then, with a constant Co = 
C (a Y ,T J), 

sup \\u e - w £ \\ L 2 {R n) +La o(^ n) < C e . (1.6) 
te[o,T e- 2 } 

The definition of the norm in (1.6) is recalled at the end of Section 2.2. We 
note that / G L 1 (IR n ) implies F G L°°(IR n ). Since F has compact support, we 
obtain F G L 2 (IR n ), / G L 2 (IR n ), and the regularity / G H s (R n ) for every s > 0. 
From now on, we understand the Fourier transform (1.3) in the sense of L 2 (IR n ). 

Comparison with the literature 

The derivation of effective equations in periodic homogenization problems is an 
old subject [22], two-scale convergence [1] is today the most relevant analytical 
tool. The use of Bloch-wave expansions was explored only more recently, see 
e.g. [7]. 

Compared to elliptic and parabolic equations, some distinctive features are 
relevant in the analysis of the wave equation. One observation of [4] was that 
convergence of energies can only be expected for initial data that are adapted to 
the periodic medium, see also [14]. Diffraction and dispersion effects are analyzed 
in the spirit of homogenization theory in [2, 3]. While the underlying questions 
are similar, these contributions study a different scaling behavior in e. Other 
homogenization results for the wave equation are contained in [5, 15, 19, 20, 24, 25]. 

The interest in the derivation of a dispersive effective wave equation appears 
most clearly in the works of Chen, Fish, and Nagai, e.g. [10, 11, 12, 13]. The 
authors expand several ideas to treat the problem, among others they propose to 
introduce a slow and a fast time scale to capture the long-time behavior of waves. 
The authors concentrate on numerical studies and do not provide a derivation of 
an effective model. 

Derivation of dispersive models. To our knowledge, the first rigorous result 
that establishes a dispersive model for the wave equation in the scaling of (1.1) 
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appeared in [18]. In that contribution, the one- dimensional case was analyzed, 
the one-dimensional version of (1.5) was formulated (in this case, A, E, and F 
are scalar coefficients and the differential operator is D = d x ), and a result similar 
to our Theorem 1.1 was shown: the well-posedness of the dispersive equation and 
an error bound on large time intervals. 

Beyond the one dimensional case, we are not aware of any rigorous results. The 
most relevant contribution with the perspective taken here is [23]. In that paper, 
Bloch-wave expansions are used to analyze the problem, mathematical insight is 
gained, and the dispersive wave equation (3.1) is formulated (not in one of the 
theorems, but as a formal consequence on page 992). We use many of the ideas 
of that contribution. 

Equation (3.1) appears also as equation (42) in [13], the authors call it the 
"bad" Boussinesq equation. The problem about this equation is its ill-posedness: 
Loosely speaking, the equation is of the form d 2 u + Lu = 0, with L = —A — eA 2 . 
The lowest order part (in e) of L is —A, hence a positive operator, but for every 
e > 0, the operator is negative, since A 2 is positive and contains the highest order 
of differentiation. One can speculate that this was the reason why no effective 
dispersive models were rigorously formulated in the above mentioned works. 

It was already observed in [13], that a "good" Boussinesq equation can be 
obtained with a simple trick: Loosely speaking, in an equation of the form d 2 u = 
—Lu = Au + eA 2 u, we can replace Au to lowest order (in e) by d 2 u, hence we 
may write the equation as d 2 u = Au + eAd 2 u. In this form, the equation is well- 
posed. This observation was also exploited in [18], where it was shown rigorously, 
that the "good" Boussinesq equation is the effective model for large times in the 
one-dimensional case. 

In this contribution we treat the higher dimensional case. The methods are 
completely different from those of [18], but are based on Bloch-wave expansions 
as used in [23]. We must introduce two assumptions: (i) we use inital data that 
are compactly supported in Fourier space and (ii) we assume some symmetry of 
the medium. Under these assumptions, we derive an effective dispersive equation. 
Regarding the replacement that transforms the "bad" effective equation into a 
"good" one, we must work with tensors of coefficients, but the essential idea 
remains the same. We show with mathematical rigor that the new equation has 
the desired approximation property. 

In Section 3 we expand the solution u e in Bloch waves, in Section 3 we analyze 
the weakly dispersive equation (1.5). The proof of Theorem 1.1 is concluded at 
the end of Section 3. 

2 Approximation with a Bloch wave expansion 

In this section we present, in slightly changed notation and with mathematical 
rigor regarding assumptions and norms, the approximation results of [23]. To 
simplify some of the notation of [23], we consider here only the mass-density p = 1 
and the scaling factor A = 1. We abbreviate the square of certain eigenvalues with 
/i := U) 2 . 
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2.1 Bloch wave expansion 

We are given a medium by the coefficient dy(y) on the cube Y. The Bloch 
wave expansion uses functions if) m , which are solutions of a periodic eigenvalue 
problem on Y . The wave parameter k is a vector in the reciprocal periodicity 
cell Z = (— 1/2,1/2)". At this point, we regard k G Z as a given parameter and 
consider 

- (V y + \k) ■ (a Y (y)(V y + ik)ip m (y, k)) = fA m {k)ip m (y , k) . (2.1) 

We search for i/) m (.,k) : Y — > C in the space H^ cr (Y), defined as the space of 
periodic functions on Y of class H 1 . We find a family (indexed by m G N = 
{0, 1, 2, . . . }) of periodic solutions 4> m (., k) : Y — )■ C with eigenvalues fi m (k), both 
the solution and the eigenvalue depend on k. We assume that the functions are 
normalized in L 2 (Y), \\ip m \\L 2 (Y) = 1- Regarding the regularity of ip m we note 
that, for ay of class C 1 , standard elliptic regularity theory implies ip m G H 2 (Y). 

Based on the eigenfunction if) m , we can construct the quasi-periodic function 
w m{y, k) := ip m {y, k)e lk y . The eigenvalue problem (2.1) is designed such that w m 
has the eigenfunction property 

- V y • (a Y (y)V y w m (y, k)) = fi m (k)w m (y, k). (2.2) 

We recall an essential fact regarding the completeness of these eigenfunctions 
(see e.g. [7] for this well-known result). The Bloch solutions form a basis of 
L 2 (M. n ) in the sense that every function g G L 2 (IR n ) can be expanded as 



oo „ „ 

(v) = ^2 / 9m(k)w m (y,k)dk, g m (k) = \ 

m=0 Jz jRr 



g(y)w m (y,ky dy , (2.3) 



where we use the star * to denote complex conjugation. There holds the Parseval 
identity 

p OO p 

Iblli^K") = / \9(x)\ 2 dx = Y" / \g m (k)\ 2 dk= \\g\\pm L 2 {z)) ■ (2.4) 
Rescaled Bloch wave expansion 

We investigate a strongly heterogeneous medium a £ (x) = ay{x/e). Starting from 
the Bloch waves on the cube Y, we therefore define rescaled quantities as 

i>m( x > k ) '■= (~,ek) , n e m {k) ■= \nm(ek) , (2.5) 



e 2 ' 

<(x, k) := w m ek) = ^(x, k)e ik ' x = ^ m ek) e ik x . (2.6) 

This choice guarantees, in particular, 

- V • (a%x) V<(x, k)) = li E m (k)w £ m (x, k). (2.7) 
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The expansion formula (2.3) in Bloch eigenfunctions can be expressed in the 
new variables. Every function / G L 2 (M. n ) can be written as 

/(*) = E / k) dk , t(k) = / f(x)w £ m (x, k)* dx . (2.8) 

m=0 J Z/* n 

To verify this formula, it suffices to set f(x) = g(x/e) and f^,(k) = e n g m (ek). 
This shows additionally the Parseval identity in transformed variables, 



2 



.,= / \f( X )\ 2 dx = Y, \fm\ 2 dk = \\ni [W z/e))- (2-9) 



Expansion of the solution 

The Bloch-wave formalism can provide a formula for the solution. 

Lemma 2.1 (Expansion of the solution). Let the Y -periodic coefficients be of 
class ay G C 1 (M ?1 ) and positive, let initial data f G H^iW 1 ) be given as in (1.3). 
Then, for every e > 0, the wave equation (1.1) has a unique weak solution u e with 
the regularity u £ (x,t) e L°°([0, oo), H 2 (R n )) n H' rl,oo ([0, oo), H 1 ^ 71 )). 

We expand the initial values f as in (2.8) and make the assumption that the 
series is convergent in f/" 1 (R n ). Then the weak solution u e of (1.1) can be repre- 
sented as 

oo „ 

u £ (x,t) = V/ f^(k)w £ m (x,k)Re(e it '/^ ir) ) dk. (2.10) 

The right hand side is understood as the strong L 2 (W n ) limit of partial sums, for 
every fixed t > 0. Re(.) denotes the real part. 

Before we start the proof, we note that the expression in (2.10) formally defines 
a solution of (1.1)— (1.3). In fact, the second time derivative of the right hand side 
is given by the same formula, introducing only the additional factor —if m {k) under 
the integral. On the other hand, the application of the operator V • (a e (x)V) to 
the integrand produces, by (2.7), the same result. 

Proof. Step 1. The energy solution. The solution u £ can be constructed, e.g., 
with a Galerkin scheme. One exploits the energy estimate which is obtained with 
a multiplication of the equation with (the real function) dtu £ 



0=[ [d 2 u £ (.,t)-V-(a £ Vu £ (.,t))}d t u e = ~ [ \d t u £ \ 2 + a £ \Vu £ (.,t)\ 2 . 

Also higher order estimates can be obtained. We use L £ := V ■ (a £ (x)'V) and 
multiply the equation d 2 u £ = L £ u £ with —dt(L £ u £ ) to find 

f t \ f a £ \d t Vu £ \ 2 + \L £ u £ \ 2 = 0. (2.11) 

Since the initial data are u(t = 0) = / G i? 2 (lR n ) and dtu(t = 0) = 0, we obtain 
estimates for u £ in the stated function spaces. Estimates for L £ u £ G L 2 imply 
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estimates for u e G H 2 due to ay G C 1 by standard elliptic regularity theory. 
Uniqueness within the given class follows from linearity and a similar calculation. 

Step 2. Convergence in (2.10). The Parseval identity (2.9) implies that the 
coefficient functions define an element (fm(^))m,k of / 2 (N, L 2 (Z/e)). As a conse- 
quence, also the modified coefficients (/^(fc)Re ( e^V^w) ) define an element 



m,k 

in the same space, since all factors have the unit norm. Using again the Parseval 
identity (2.9), we conclude that the sum of (2.10) converges in L 2 (R n ), indepen- 
dent of t > 0. 

Step 3. Identification of u £ . We consider a partial sum Ylm=i ^ n (2- 10) to 
define a function u e M and observe that this provides a strong solution u £ M of the 
wave equation to the initial values /m = J2m=o Iz/e fm(k) w m( x i ^) dk and van- 
ishing initial velocity. This fact can be checked with the above-mentioned formal 
calculation, the operator V ■ (a e (x)V) is understood in the weak form and can 
be applied to the if 1 (F)-functions w m . We claim that u £ M forms a Cauchy se- 
quence in every space L°°([0, T], i/ 1 (IR n )). This follows with a testing argument, 
exploiting 

/ a £ \Vul I (t)-Vu%(t)\ 2 +\d t ul I (t)-d t u%(t)\ 2 = [ a £ \Vf e M - V/^| 2 ^ 

for M,N — > oo due to the assumed if 1 — convergence of (2.8). We conclude that 
u e M converges to a limit function. The limit function is again an energy solution 
of the wave equation, by uniqueness we conclude u e M — > u e for M — > oo. 

On the other hand, as observed in Step 2, by definition of u £ M , the limit function 
is given by the right hand side of (2.10). □ 



2.2 The approximation results of Santosa and Symes 

With the next two theorems we observe that, for small e > 0, the expression of 
(2.10) may be simplified. In our first simplification we realize that all indices m 
with m > 1 can be neglected. 

Theorem 2.1 (Santosa and Symes [23], Theorem 1). Let ay be of class C 1 , let 
the initial values be given by f G H 2 (M. n ) of the form (1.3), such that the series 
of (2.8) is convergent in ^(W 1 ). Let u £ G L°°((0, oo); H 2 (W)) be the solution of 
(1.1), given by (2.10). Then there exists C = C(f) > such that 



sup 

tg(0,oo) 



oo „ 

E/ L(k)w* m (x,k)ReC 

m=l J Z/e V 



3 i*\/Vm( fc ) 



dk 



< Ce 



(2.12) 



L 2 (R") 



Proof. We consider a single coefficient /^(/c)e ±1 * in the expansion of u £ in 
(2.10). We use first the inversion formula (2.8) to express this coefficient, then 
the eigenvalue property (2.7) to introduce the factor /^^(fc), then integration by 
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parts and the solution property of u £ , 

fa(k) Re (eV/^T) = ! u £ (x,t)w £ m (x,k)*dx 

= — TTJa I u £ (x,t)[V ■(a £ (x)Vw £ m (x,k))]*dx 

= JTiaI ld?u £ (x,t)]w £ m (x,kydx. (2.13) 

We claim that, with C > independent of t £ [0,oo), the functions x h> 
d 2 u £ (x,t) satisfy the estimate ||<9 2, u e (., £)||z,2( R n) < Ce~ x . Indeed, this bound can 
be obtained as in (2.11), where multiplication of d 2 u £ = L £ u £ with dtL £ u £ provided 

/ a £ \d t Vu £ (.,t)\ 2 + \L £ u £ (.,t)\ 2 = [ |LV(,0)| 2 . 

Since initial data / are smooth, we have ||ij e « £ |t=o||.L 2 = || V-(a e (x)V/)||L2 < Ce~ x , 
hence \\L £ u £ (., t) \\ L 2 < Ce~ x . Accordingly, by the evolution equation, we also have 
\\d?u £ (.,t)\\ L2 = \\L £ u £ (.,t)\\ L2 < Ce- 1 . 

We can now continue (2.13). Using again the Parseval identity, we conclude 



k)f m (k)Ke (e^V^m) = e 2 \\d 2 u £ (.,t)\\ L2{Rn) < Ce . 

\ / i^{ri,L z (Z/e)) 



It remains to observe that omitting the term m = decreases the norm on the 
left hand side. Regarding terms with m > 1, we exploit that there exists a lower 
bound Co > such that /x m (0 > Co, independent of £ £ Z and m > 1, cf. [7]. 
Another application of the Parseval identity provides the claim (2.12). □ 

At this point, we have obtained a first approximation of the solution. In the 
expansion of u £ , all contributions from indices m > 1 are not relevant at the lowest 
order (uniformly in time). Theorem 2.1 thus reads \\u £ — ^olk 00 ((o,oo),L 2 (R n )) — Ce, 
where 

u £ (x,t):=J f £ (k)w £ (x,k)Re(e u ^^^ dk . (2.14) 

Our next aim is to replace the Bloch coefficient f${k) by the Fourier coeffi- 
cient Fq. At this point, we make more substantial changes with respect to [23]. 
Nevertheless, the essence of the subsequent result has been observed in Theorem 
2 by Santosa and Symes, [23]. 

Theorem 2.2. Let the dimension be n > 1, let the initial data f £ L 2 (R") n 
L 1 (M n ) satisfy (1.3) with F £ L 2 (IR n ) supported on the compact set K. We 
assume that the coefficient ay has a regularity that implies ipo(-,Q £ C l (Y,M), 
with a bound that is independent of £ £ Z. Then, with C = C(f) > 0, there holds 



<Ce. (2.15) 

LHZ/e) 



fo-Fo 

Furthermore, fore > small enough to satisfy diam(J^) < e" 1 , there holds 

f^k)=0 Vk£(Z/e)\K. (2.16) 
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We note that the regularity requirement on ip is satisfied, e.g., for ay G C 1 , 
see the proof of the main theorem at the end of Section 3. 

Proof. Step 1: k G K . We write the difference of the two functions for fixed k as 

1 



m)-F (k) 



-ik-x 



[ ^,ek 



dx . 



The periodic solution ?po(.,0) to the wave vector k = is constant, by our nor- 
malization it is given as ijjo(y,0) = \/\Y\ for every y G Y. Since k ranges (in 
this step of the proof) in the bounded compact set K, we find the estimate 



< Ce 



(2.17) 



uniformly in k G K and x G W 1 for some constant C = C(ay). This is a 
consequence of the fact that ^ (->0 £ C°(Y) depends in a Lipschitz continuous 
way on £. We assumed / G L 1 (R n ) and obtain therefore 



m)-F (k) <Ce\\f\\ LH 



< Ce 



uniformly in k G K. Since K is compact, this provides also an L 1 (i^)-bound as 
in the statement of (2.15). 

Step 2: k G Z/e \ K. We expand ^o(-;0 m a Fourier series on the cube Y, 
treating £ clS db parameter, 

lfo(l/,0 = £>(0e U " 1 '- 
Since we assumed that ipo(-, is of class C 1 , the Fourier series converges uniformly 



in y, for every parameter £. Since / is of class L 



we may write 



f(x)e iIjq ( ~,ek) dx 



x 



lim V oi{ek) [ f(x)e-' lkx e iHx/£) dx. 



iez n ,\i\<L 

For each single term we find, since k — (l/e) G" K for Z n 9 Z 7^ 0, 



/(x)e- ifc -V' (x/£) dx = (27r) n/2 F (A; - (Z/e)) 







(2vr)™/ 2 F (fc) 



for Z 7^ 
for Z = 



Evaluation of the sum, we find 

/ £ (fc)=«o(^) (2vrr/ 2 F (fc) = 0. 

In particular, we obtain the claim (2.16) about the support of /q. This, in turn, 
implies also the L 1 -estimate (2.15) for the difference on all of Z/e. □ 
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We want to use the theorem in order to simplify the representation of the 
solution Uq of (2.14). We define a new approximation as 



U £ (x,t) := (2ny n/2 



J Fo^e^Re^V^) dk. (2.18) 



Theorem 2.2 allows to calculate, using once more (2.17) to compare Wq(x, k) 
ifjo(x/e,ek)e ik - x with (2vr)- n / 2 e ifc - x , 



\ U ~ ^ £ ||l°°((0,oo)xR™) 



K 



f £ (k)w £ (x, k) Re ( e i4 vW) ) dk - U 



L°°((0,oo)xR n ) 



< 



-\— sup sup / f^k^Refe^^m) dk 
-{2-K)- n/2 J F (k)e ik x Re (V'v 7 ^)) dk 



< C 



fo-Fo 



+ Ce< Ce. 

LHZ/e) 

We can combine this error estimate with the one obtained earlier for the difference 
\\ u£ — n olU°°((o,oo),L 2 (iR"))- We use, given two norms ||.||x and ||.||y, the new norm 
(weaker than both original norms) ||it||x+y := inf{||ui||.x + ll u 2||y : u = u± + u 2 }. 
This allows to write the combined estimate as \\u e — ^ £ ||l° o ((o,oo),(l°°+l 2 )(K")) < Ce. 



2.3 Expansion of the dispersion relation 

The next step is to replace the eigenvalue /j,q by its Taylor series. We note that in 
a neighborhood of k = the eigenvalue /io depends analytically on k and /Uo(0) = 
V/io(0) = 0, cf. [7]. We denote the derivatives of fi via A im = ^d kl d km fi (0), 
Bi mn = ld kl d km d kn n (0), and C lmnq = ^d kl d km d kn d kq n (0). We will below study 
a symmetric situation in which B vanishes. We can then assume that the Taylor 
series of /io in k around k = is given as 



Ik) 



(k) = A imkik m + C lmnq kik m k n k q + 0(\k\ 5 ). (2.19) 



Here and below a bare sum is always over the repeated indices. The expansion 
corresponds to the following expansion of fJ%(k), 

(*0 = ^o( £k ) = ^2 A im k ik m + e 2 Ci mnq kik m k n k q + 0{e 3 ) , (2.20) 

the error is of order e 3 , uniformly in k G K. 

In the spirit of this expansion, we next want to simplify further U e of (2.18). 
We define v e (compare page 992 of [23]) as 



v E {x,t) := (2n)- n ' 2 



\Y.j K F o(k)e ik - x exp (iity^A^fc^ 



x exp ±— t 



, C\ mnqkik m k n k q 



(2.21) 



2 y/J2 Aimkk 



dk 
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We arrive at the following approximation result. We repeat that the underlying 
observations are taken from [23]. 

Corollary 2.1. Let f G L 1 (R n ) have a Fourier-transform F supported on a 
compact set K C M n . Let u £ be the solution of (1.1) and let v £ be defined by 
(2.21). Then 

SUP \\u £ (t) - ^(t)||L2 (R n ) + i oo (K n } < Cs. (2.22) 

t£{0,T e~ 2 } 

Proof. The estimate for the difference u £ — U e has been concluded after the defi- 
nition of U £ in (2.18). It remains to estimate the difference v £ — U e in the same 
norm. 

The Taylor expansion of the square root reads 

Va + c = Va + -^-=c + 0{\c\ 2 ). 
2y a 

This implies that the definitions of U £ (t) and v £ (t) coincide, except for a factor of 
the form 

exp (±itO(e 4 )) = l + 0(e 2 ), 

uniformly in t for t G [0,To£ _2 ]. Because of Fq G L°°(IR n ) and boundedness of K, 
this implies (2.22). □ 

In view of Corollary 2.1, it will no longer be necessary to work with u £ , the so- 
lution to the original wave equation in a heterogeneous medium. We will, instead, 
restrict ourselves to the analysis of the function v £ defined by (2.21). 

Note that Taylor expansions of Bloch eigenvalues are commonly used also in 
the derivation of effective equations for envelopes of nonlinear waves in periodic 
structures, see e.g. [8, 9]. 



2.4 Symmetries 

The structure of the three tensors A, B and C, defined via the expansion of Ho{k), 
is very simple if we consider only symmetric material functions ay- Indeed, we 
will see that A, B, and C are fully characterized by three real numbers a*, a, and 
P. 

We assume that ay is symmetric with respect to reflections across a hyperplane 
{Uj — 0}, j G {0, ...,n}, and invariant under coordinate permutations. To be 
more precise, we introduce the following transformation of lR n , defined for y = 
(y h ...,y n ) as 

Si(y) = (Vi, ■ ■ ■ , Vi-i, -Vi, Vi+i, ■■■,Vn), 

R ij(y) = (yi, ■ ■ ■ > i/i-i, yj, ■ ■ ■ , Vj-u vu Vj+i, •■■,!/«)■ 

Our symmetry assumption on ay can now be formulated as 
a Y (y) = a Y (Si{y)) = a Y {Rij{y)) for all i,j G {1, . . . ,n} and all y G M n . (2.23) 
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As we show next, the symmetry properties of ay in y imply the identical symmetry 
properties of /io in k, 

/x (fc) = ^o(Si(k)) = fi (Rij(k)) for all i, j G {1, . . . ,n} and all k G K. (2.24) 

In fact, (2.24) holds also for all functions fi m , but we exploit here only the symme- 
try of fio- To show (2.24), we express fio(k) with the variational characterization, 
see Theorem XIII. 2 in [21], as 

fio(k) = min I(w,k), where I(w, k) := / ay(y)|(V + ik)w\ 2 dy . (2.25) 

™e-ffper<T) Jy 

Using the symmetry of ay, we can calculate 

I(w, Si(k)) = J^a Y {y) |[(V + iS t {k)M iv)? dy 

= I a Y {y)\[{V + iS i {k))w]{S l {y))\ 2 dy (2.26) 

JS-\Y) 

= a Y (y) IS, ([(V + ik)(w o S^] (y))| 2 dy = I(w o S u k) . 

Minimizing over the functions w o Si provides the same result as minimizing over 
w, since with w G H^iY) also wo Si G H^ er (Y). This provides (2.24) for Si. The 
calculation for Rij is identical. 

As a consequence of the symmetry, we obtain the following characterization of 
the Taylor expansion coefficients. 

Lemma 2.2. Let ay satisfy the symmetries (2.23). Then the tensors A, B and 
C, defined in (2.19), satisfy 

An = A u =: a*, A {j = 0, B = 0, 

Can = Cim =: a, Cyij = Wjjj = Cujj = C1122 =: ft 

for all i,j G {1, . . . , n} with i ^ j. All entries of C , that are not mentioned above, 
vanish. 

Proof. The proof uses symmetry (2.24). The symmetry under Si implies that fio 
is an even function. Thus all derivatives of /io with an odd number of derivatives 
in one variable vanish at k — 0. This proves A^ = 0, B = 0, and, e.g., Cmj = 0. 
The fact that derivatives can be interchanged provides, e.g., Cnjj = Cmj. 
The symmetry under allows to calculate 

dlMk) = <9 fe 2 >o o R^ik) = [dl^]{R l3 {k)) . 

Evaluating in k = provides An = Ajj. The analogous calculation for fourth 
derivatives shows, e.g., Can = Cjjjj. This proves the claim in the two-dimensional 
case. 

For n > 3 we can calculate again with the symmetry under R^ 

9lAM k ) = <t^>° = [9lMM( R M)- 

Evaluating at k = 0, we obtain Cnjj = Can for all indices i, j, I < n. □ 
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3 A well-posed weakly dispersive equation 

A weakly dispersive equation that is related to the definition of v e is (at this point, 
we correct a typo of [23] regarding the sign before C) 

d 2 u = AD 2 u-e 2 CD 4 u. (3.1) 

Indeed, when applied to v e , the operator AD 2 produces the factor —Ai m kik m under 
the integral, and the operator —e 2 CD 4 produces the factor —e 2 Ci mnq kik m k n k q . 
The second time derivative produces the factor 

A\ m k\k m £ Ci mn qkik m k n kq (e /^){Ci mn qkik m k n k^) /(Aimkikm) 

under the integral. Therefore, up to an error of order e 4 , the function v e solves 
(3.1). 

We emphasize that, in general, (3.1) cannot be used as an effective dispersive 
model. The fourth order operator —CD 4 on the right hand side can be positive 
such that (3.1) is ill-posed. In the one- dimensional setting, C < is shown in [18], 
hence the equation is necessarily ill-posed. Section 4.2 includes a two-dimensional 
numerical example where the numbers a and /3, describing C, satisfy a < and 
> 0. Moreover, there holds 3/3 < \a\, such that —CD 4 is a positive operator. 

As a consequence, even though v e solves (3.1) up to an error of order e 4 , we 
cannot conclude that solutions to this equation provide approximations of v e . 
Even worse, it may be impossible to construct any solution of (3.1). 

3.1 Decomposition of the operator for symmetric media 

As indicated in the introduction, our aim is now to replace (3.1) by a well-posed 
equation, which is equivalent in all relevant powers of e. We therefore start from 
the two tensors A = a* id G R nxn and C G ]& nxnxnxn f Lemma 2.2 and consider 
the operator 

n n 

CD 4 = ]T C l3k Ad 3 d k d { = a df + 3/3 ]T d 2 d 2 . (3.2) 

ijkl i=l «>j=i 

To avoid confusion, we note that = ^J2i<j- O ur a i m i s to construct coef- 

ficients E G M nxn and F G ]R nxnxnxn suc fi that the differential operator can be 
re-written as 

— CD 4 = ED 2 AD 2 — FD 4 , (3.3) 
where E and F are positive semidefinite and symmetric, i.e. 

n 

F W&jtki > for every £ G M nxn and Fijkl FfaUj (3.4) 

i ,j,k,l=l 

and Y^j=i EijViVj ^ f° r every 77 G M. n and = Eji for i,j,k,l G {l,...,n}. 
The decomposition result (3.3) allows, using the lowest order of (3.1), to re- write 
the operator in the evolution equation formally as 

- e 2 CD 4 u = e 2 ED 2 AD 2 u - e 2 FD 4 u = e 2 ED 2 d 2 u - e 2 FD 4 u + 0{e 4 ) . (3.5) 

With this replacement in equation (3.1), we obtain the well-posed equation (1.5). 
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Lemma 3.1 (Decomposability). Let A £ R nxn and C £ ]R« X ™ X ™ X ™ be as in 
Lemma 2.2, given by three constants a* > 0, a,/3 £ R, in particular with CD 4 
given by (3.2). Then there exist symmetric and positive semidefinite tensors E £ 
R nxn and F £ R™x™ x ™><™ ^ ^ 4 ^ ^ wntten as m (3.3). 

[/sing {a} + := max{a, 0} to denote the positive part of a number a, a possible 
choice of E and F is 

E u = \ {{-<*}+ + 3{-P}+), E tJ =0, (3.6) 

(X 

Fan = {a} + + 3{-f3} + , F m = {-«}+ + 3{/3} + , (3.7) 
for all i,j £ {1, ...,n} with i 7^ j. All other entries of F are set to zero. 
With (3.6)-(3.7), we introduce the two differential operators 

n 

ED 2 = - ({-a} + + 3{-/3} + ) J2% = - + 3 {-^}+) A ' 

a* £ r^ a* 

1=1 

n n 

FD 4 = ({«} + + 3{-/3} + )^^ + ({-«} + + 3{/3} + ) d 2 d). 

i=l 

Since a and (3 are real numbers, there are four different possibilities for the 
signs of a and /?. Distinguishing these four cases, we can write the two differential 
operators in very simple expressions. 

Remark 3.1. The operators ED 2 and FD 4 of (3.6)-(3.7) are given as follows. 
Case 1. a < 0, /3 < 0: 

n n 

ED 2 = — (lal + 3|/3|)A and FD 4 = 31/31 V df + \a\ V d 2 d 2 
a* ^— ' ^— ' J 

i=i ij=\^j 

Case 2. a < 0,/3 > 0: 

1 1 n 
-E-D 2 = and FD 4 = (\a\ + 3/3) ^ <9 2 <9 2 . 

Case 3. a > 0,(3 < 0: 
ED 2 

Case 4. a > 0, /3 > O.- 
ED 2 = and FD 4 = a^9f + 3/3 d 2 d 2 = CD 4 . 



U and FD 4 = (a + 3|/3|)^<9 4 
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We note that the first two cases ( with a < 0) are the relevant ones in our numerical 
examples. 

Proof of Lemma 3.1. Step 1. Properties of E and F . By definition, E is a nonneg- 
ative multiple of the identity in M n . The tensor is therefore positive semi definite 
and symmetric. Also F is symmetric by definition. For £ e M. nxn holds 

n 

/] Fijki&jCki 

i,j,k,l=l 

n n 

= E(M + + 3{-/?} + )fe) 2 + Yl ({-«} + + 3{/3} + )te,) 2 >0. 

Hence F is also positive semidefinite. 

Step 2. Decomposition property. It remains to show —CD 4 = ED 2 AD 2 —FD 4 . 
For that purpose we calculate the right hand side 

ED 2 AD 2 - FD 4 

n / n \ n 

= - a-a} + + 3{-/?} + ) E 9 n E a * d * - (M+ + E d t 

1=1 \j=l / 1=1 



-({-«}+ + 3{/3} + ) £ 



2 



({-«}+ + 3{-/?} + )J>f + ({-<*}+ + 3{-/3} + ) ^ 



J 2 

i=i ij=i t i^j 

n n 

-(W + + 3{-/3} + )^^-({-«} + + 3{/3} + ) £ «9 2 «9 2 
i=i ij=i,i^- 

n n 

i=i 1^=1^ 

This is the desired decomposition (3.3). □ 



3.2 An approximation result 

With the subsequent theorem, we provide the central error estimate for our main 
result. We start from two tensors A and C (in the application of the theorem they 
are defined by (2.19)), and assume that C is decomposable with tensors E and F. 
With these four tensors we can study two objects: The solution w e of (1.5), and 
the function v e , defined by the representation formula (2.21). Our next theorem 
compares these two objects. 

Theorem 3.1. Let A,C,E,F be tensors with the properties: A is symmetric 
and positive definite, -^ijCiCj — 7l£| 2 f or some 7 > 0, E and F are positive 
semidefinite and symmetric, C allows the decomposition (3.3). Then the following 
holds. 
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1. Well-posedness. Let R G L 1 (0, T e~ 2 ; L 2 (M n )) be a right hand side and 
let f G if 2 (IR n ) be an initial datum. We study an inhomogeneous version of 
equation (1.5), 

d 2 w £ (x, t) - AD 2 w £ (x, t) - e 2 d 2 ED 2 w £ {x, t) + e 2 FDV(x, t) = R(x, t) , 
w £ (x,0) = f(x), d t w £ (x,0) = 0. 

(3.8) 

for x G M n and t G (0,T £~ 2 ). This equation has a unique solution w £ G 
L co (0,T e~ 2 ;H 2 (W 1 )). 

2. Approximation. Letv 6 be defined by (2.21) with Fq and f related by (1.3). 
Let w £ be a solution of (3.8) to R = 0. Then 

SUP \\d t (v e - W £ )(.,t)|| L 2 (M n ) + SUP ||V(l> £ -W £ )(.,t)|| L 2 (R n ) < C £ 2 , 

tG[0,T e- 2 ] te[0,T () e- 2 } 

(3.9) 

where C > denotes a constant that depends on f and the coefficients, but 
is independent of e. 

Proof. Well-posedness of problem (3.8). We use the following concept of weak 
solutions. We say that w £ G L°°(0, T e~ 2 ; H 2 (W n )) with the property d t w £ G 
L°°(0, T Q e~ 2 ; if 1 (M n )) is a weak solution, if it satisfies w £ (x, 0) = f(x) in the sense 
of traces and if 



pT e 2 p pT e - p 

/ R(p= / {-dtw e dt<f> + V0- AV#} 



5 



rT Q e- 2 p 

1 / / {-V(<9 4 0) • £V(3 t w e ) + D 2 : FD 2 w £ } 

JO JR n 



(3.10) 



for every test-function <p G C£([0, T £- 2 ); # 2 (IR n )). Here L> 2 : FD 2 w e denotes 
the tensor product of D 2 and FD 2 w £ , 

D 2 ct>:FD 2 w £ := did^F^dkdtw'. 

i,j,k,l=l 

We prove the existence of a weak solution to problem (3.8) with a Galerkin 
scheme. We use a countable basis {ip k }km of the separable space H 1 (M. n ) and 
the finite-dimensional sub-spaces Vk '■= spanfy 1 , ...,ip K } C // 1 (lR n ). The basis 
{i/j k }keN is chosen in such a way that the functions ip k are of class H 2 (M. n ) and 
such that the family of L 2 -orthogonal projections Pk onto Vk are bounded as 
maps Pk '■ i/ 2 (M n ) — > H 2 (M. n ). For every K G N we search for approximative 
solutions w;^ of the form 



w K 



K 

£ 



[0,T e- 2 ]-*V K , = 5>*W 



fc=i 
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with coefficients a| : [0, T$e 2 ] — > M.. We demand that w e K solves (3.8) in the weak 
sense, however, only for test-functions in the i^-dimensional space Vk, 

R^ k = f {d 2 t w £ K ^ k + V^ fc • AV^} 

Jr " (3.11) 
+ e 2 / {Vi) k ■ EV{d 2 w e K ) + D 2 ip k : FD 2 w £ K } 



for every k G {1, K}. For the initial data we demand that (w e K \ t =o, ip k ) L 2 (R n ) = 
(/, ^ fc )L 2 (M") and (d t w £ K \ t =o, ip k )^(R") = 0. For every K e N, equation (3.11) is 
a i^-dimensional system of ordinary differential equations of second order for the 
coefficient vector (a\ (t), . . . , a e K (t)) , which can be solved uniquely. This provides 
the approximative solutions w £ K . 

We now derive i^-independent a priori estimates for the sequence w e K . For 
that purpose we test equation (3.8) with d t w e K (more precisely, we multiply (3.11) 
with d t a e k and take the sum over k). Exploiting the symmetry of A,E and F we 
obtain 



/ Rd t w K =\d t [ {\d t w £ K \ 2 + Vw K -AVw K } 

Jr" * JR n 

+ z 2 \dt [ {V(<9 t <) • EV(d t w K ) + D 2 w s K : FD 2 w e K }. 



(3.12) 



We next integrate (3.12) over [0,t Q ], where to £ [0,To£~ 2 ] is arbitrary. We exploit 
the initial condition w e K \ t= o = fx, where fx is the L 2 -projection of / onto Vk- 
The other initial condition is d t w K \ t= o = and we arrive at 

to 

R d t w £ K + I Vf K -AVf K + e 2 I D 2 f K :FD 2 f K 



n 



= / {\d t w K \ t=t J 2 + Vw K \ t=t0 ■ AW K \ t=t0 } 

JR n 

+ e 2 [ {V{d t w K )\ t=t0 ■ EV{d t w K )\ t=t0 + D 2 w K \ t=t0 : FD 2 w K \ t=to } 

> ll^x(^o)|li 2 (Rn ) +7||V^(.,to)|li 2(R .). (3.13) 

In the last line we exploited that A is positive definite with parameter 7 > and 
that E and F are positive semi-definite. Introducing Y(t) := \\d t w K {., t)\\ 2 L2( ^ n ^ + 
7l|Vw4(., ^H^n) for the right hand side of (3.13) and Y := J* Rn {V f K ■ AVf K + 
e 2 D 2 fx '■ FD 2 f K }, we can calculate with the Cauchy-Schwarz inequality 



Y(t)<2 [ \\R(.,s)\\ L 2 {Rn) \\d t w K (.,s)\\ L 2 {Rn) ds + Y 
Jo 

<2 / \\R(.,s)\\ L 2 { Rn)VY^jds + Y . 
Jo 



(3.14) 



We claim that a Gronwall-type argument leads from inequality (3.14) to the esti- 
mate 

Y(t)<2Y + 2( [ \\R(.,s)\\ L 2 {Rn) ds) . (3.15) 



18 



T. Dohnal, A. Lamacz, B. Schweizer 



For the sake of clarity we postpone the justification of this implication to the end 
of the proof. With inequality (3.15) at hand we finally obtain the following a 
priori estimate 

sup Y(t) = sup {||^(.,t)||i 2(R „ ) +7||V^(.,t)||i 2(R „ ) } 
te[o,T 0£ - 2 ] te[0,T o e- 2 } L } 

< 2Y + 2||i?|||i( 0)To£ - 2 . i 2( R „)) (3.16) 

< 2(C(A) +e 2 C(F))\\f\\ 2 H2{Rn) + 2\\R\\ 2 Ll{0tTo£ ^ L2(my 

The bound in (3.16) is independent of K. Hence, possibly after passing to a 
subsequnce, we may consider the weak limit K — > oo of solutions w £ K of the 
Galerkin scheme. Due to the linearity of the problem, the limit provides a solution 
w £ E L°°(0,T £- 2 ; H\R n )) with d t w £ e L°°(0,T £- 2 ; L 2 (W 1 )) to (3.8) in the sense 
of distributions. Furthermore, w £ satisfies exactly the same a priori estimates as 
its approximations w £ K . By differentiating (3.8) with respect to x, one discovers 
that w e has in fact higher spatial regularity and that the distributional solution 
w £ is in fact a weak solution in the sense of (3.10). Note that the uniqueness of 
solutions to problem (3.8) is a direct consequence of the a priori estimate (3.16). 
Hence, the weakly dispersive problem is well-posed. 

Proof of the approximation result (3.9). By applying the differential operator 
d\ — AD 2 — e 2 d 2 ED 2 + e 2 FD 4 to v £ , which is explicitly given in (2.21), one 
immediately discovers that v £ solves Equation (3.8) with a right hand side of 
order e A . More precisely, we calculate first with the decomposition of the operator 
-CD 4 = ED 2 AD 2 - FD A 

d 2 v £ - AD 2 v £ = -e 2 CD 4 v £ + e 4 R £ = e 2 ED 2 AD 2 v £ - e 2 FD% £ + e 4 R £ , 

where the error term comes from the double differentiation of the last factor of v 6 
with respect to time, 

R e •— _I(27r)~ n / 2 / — Clmnqklkmknkg) 2 p Q ^\ 

8 JkGK z~2 Aimkik m 

x exp | ik ■ x ± i\/y^ A m fc;fcZt | exp | ^ z [^f ^^' lmn i^ l ^ rn \ ^ 

With this preparation we can now evaluate the application of the full differential 
operator as 

d 2 v £ - AD 2 v £ + e 2 FD% £ - e 2 d 2 ED 2 v £ 
= e 2 ED 2 AD 2 v £ + e A R £ - e 2 d 2 ED 2 v £ 

(3 17) 

= e 2 ED 2 (AD 2 v £ - d 2 v £ ) + e A R £ V ' ' 

= e A ED 2 {CD\f - e 2 R £ ) + e*R £ =: R £ . 

In particular, sup tG [ OTo£ -2] ||-R £ (-, ^)|U 2 (k») < Ce A for some e-independent constant 
C . Due to the linearity of the problem and the fact that w £ is a solution to (3.8) 
with R = 0, the difference v £ — w £ solves equation (3.17) 

d 2 (v £ - w £ ) - AD 2 (v £ - w £ ) + e 2 FD 4 (v £ - w £ ) - e 2 d 2 ED 2 (v £ - w £ ) = R £ , 
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with vanishing initial data (v e — w e )(., 0) = dt(v e — w 6 )(., 0) = 0. 

By applying the a priori estimate (3.16) to the difference (v e — w £ ) we obtain 

sup |||^(^- M ; e )(.,t)||i 2(M „ ) +7||V(^-u ; £ )(.,t)||i 2(R?l) ) 
te[o,T £- 2 ] k J 

< 2 H- R£ |lLl(0,T £- 2 ;L2(]Rn)) < 2(T £ 2 ||-R e |U=o(0,Toe- 2 ;L 2 (R™))) 2 < Ce*, 

where in the last step we exploited that ||-R £ ||l°°(o r £- 2 -L 2 (M™)) < Ce 4 . This implies 
(3.9). 

Proof of the Gronw all-type Inequality (3.15). Let Y : [0, T] — > [0, oo) be a 
function such that, for a constant Y > 0, the relation 

Y{t)<2[ \\R{.,s)\\ L 2 {Rn) ^/Y^)ds + Y (3.18) 

holds for all times £ G [0, T]. We claim that then 

Y(t)<2^\\R(.,s)\\ L 2 {Rn) ds^ +2Y (3.19) 

holds for all times t e [0, T]. 

For the proof we define Z(t) to be the integral on the right hand side of (3.18), 



Z(t) :=2 I \\R(.,s)\\ L 2 m y^Y{s)ds. 
Jo 



Then Z(0) = and, due to the assumption (3.18), 

jZ{t) = 2\\R(.,t)\\ L 2 {Rn) y^}< 2\\R(.,t)\\ L 2 0ln) y/Z(t)+Y . 

We conclude that 

j t (y/Z(t) + Y ) = (2 V /Z(t) + F )~ 1 |z(t) < \\R(.,t)\\L* m . 
Integrating this relation over [0,t] we obtain, recalling Z(0) = 0, 

y/Z(t) + Y -^Y < [ \\R(.,s)\\ LHRn) ds. 

Jo 

By evaluating the square we find 

Z{t) + Y < (VYo + J Q \\R(.,s)\\ L 2 m ds 



<2F + 2 / \\R{.,s)\\ LHRn) ds 



o 



and therefore the claimed result (3.19), since Y(t) < Z(t) + Yq holds by assump- 
tion. □ 
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The main theorem. Theorem 1.1 is a consequence of the above results. 

Proof. We have seen in Lemma 2.1, that the solution u £ permits the expansion 
(2.10) in Bloch-waves. The assumptions a Y G C^Y) and / G if x (R n ) as in (1.3) 
are satisfied in the situation of Theorem 1.1. 

In Theorem 2.1 we have seen that only the term m = in the sum must 
be considered. Again, the assumptions are satisfied, ay G C 1 (F) and if 1 (!R n )- 
convergence of the Bloch-series of /. 

The assumptions of Theorem 1.1 imply / e L 2 (IR n )nL 1 (IR n ) and F G L°°(M n ), 
hence Theorem 2.2 can be applied. Furthermore, we assumed that ipo(-,£) G 
C 1 (F, R) is bounded independent of £ G This is satisfied by elliptic regularity 
theory: Holder-continuous coefficients yield bounds for solutions of divergence- 
form operators in C 1,a , see Giaquinta, [16] Chapter III, Theorem 3.2. We recall 
that Z = (-1/2, 1/2) n is bounded. 

We concluded with (2.22) a smallness result, \\u £ — v £ \\ is of order e 2 . The 
norms coincide with the ones in the claimed result (1.6) for \\u e — w £ \\, where we 
only claim the order e for the error. 

Finally, Theorem 3.1 provides the well-posedness claim and the smallness (3.9) 
of the error \\v £ — w £ \\ of order e 2 . We note that (3.9) controls only the norms of 
derivatives, but the subsequent Lemma 3.2 provides the estimate for ||t> £ — w £ \\ of 
order e in the desired norm of (1.6). □ 

Lemma 3.2. For n > 1 and T > fixed, let g £ : R n x [0,T/e 2 ] -> R be a 
sequence of functions with g £ (.,0) = 0. Then, with a constant C > 0, there holds 
an e -independent estimate 



t&[0,T/e 2 



SU P _ ||fi ,£ (-^)IU 2 (IR«)+L°°(IR") 

(3.20) 

<Ce~ l sup {\\d t g £ (.,t)\\ L 2 iRn) + \\Vg £ (.,t)\\ L 2 {Rn) } . 

te[0,T/e 2 ] 

Proof. We first consider n > 2. Given e > 0, we choose a tiling of the space as 
R n =\jE £ mJ E £ m = x m + [0,e- 1 ) n , x m = me~ 1 . (3.21) 

meZ" 

Given the function g £ we define a piecewise constant function through an averaging 
procedure, 

g £ (x,t):=-f g £ (U)dt ifxeE £ m . (3.22) 



The Poincare inequality for functions with vanishing average allows to estimate 

||<f(, t) - g%, t) = \\g%, t) - g e (,t) \\ 2 L2(E ^ 

in 

<Cdiam(^) 2 ^||V/(.,t)||i 2( ^ ) <C' £ - 2 ||V/(.,t)||l 2(M „ ) . 

m 

This provides estimate (3.20) for the part g £ — g £ . 
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In order to estimate g e , we use the fact that averaging does not increase the 
L 2 -norm, 

Y,\E e m \m{x m ,t)\^ \\d t g%.,t)\\h m < \\dtg £ (;t)\\ 2 L2{Rn) . 

With the fundamental theorem of calculus we find 

T 2 

\\9 £ (-,t)\\ 2 L°°(m.") = max |/(x m ,t)| 2 < max— sup \d t g £ (x m , s)\ 2 
m m £ se[o,T £ -2] 

^^I^J- 1 sup ^|i^||% £ 0w)| 2 
<TV" 4 sup \\d tg c(.,s)\\i 2{Rn) . 

se[0,Te- 2 ] 

For n > 2, this provides estimate (3.20) for the remaining part g £ . 

In the case n = 1 we proceed in a similar way, using now a tiling with pieces 
of larger diameter, 

R=\jE £ m , E e rn = x m + [0, e- 2 ) , x m = me-' 2 . (3.23) 

The estimate for g £ e L°°(0, Te~ 2 ; L°°(M n )) is obtained as above with the e-factor 
e~ A \E e a \~ 1 = e~ 2 as desired. To estimate the difference g e — g e we use, in the case 
n — 1, the same L°°-based norm. We calculate, for arbitrary t G (0,Te~ 2 ), 

\\g E (-:t) - <f (.,£)||l~(R") = sup||/(.,t) -^ e (.,t)IU-(^) 

■m 

< sup H^C^Uli^) < sup diam(£^) 1/2 \\d x g% t)\\^ {E ^) . 

m m 

Because of diam(i^) 1 / 2 = e -1 , this shows (3.20). We emphasize that we obtain 
a pure L°°-bound on the left hand side of (3.20) in the case n — 1. □ 



4 Numerical results 

In order to illustrate the approximation result of Theorem 1.1, we numerically 
solve equations (1.1) and (1.5) in dimensions n = 1 and n = 2 with the initial 
conditions in (1.2). 

For the spatial discretization of (1.1) we choose the fourth order finite differ- 
ence scheme of [6]. In one dimension (n = 1) and for smooth a e (x) the value of 
d x (a e (x)d x u) at the grid point x = Xj is approximated by 



22 



T. Dohnal, A. Lamacz, B. Schweizer 



where the coefficients a 6 - and a 6 , t are defined via a% = tt4— [ Xj+1 a £ (x) dx and 
= Ja/ +1 a£ ( x ) dx, and where Ax is the spacing of the uniform grid (xj)j. 
For the time discretization we use the standard centered second order scheme 
resulting in the fully discrete problem 

u m+l = 2u rn _ + ( A t) 2 (A e (A)« m )j. 

In order to initialize the scheme, we set u® = f(x 3 ) and approximate u 1 via the 

Taylor expansion u 1 = u° + A £ (\)u°. For the evaluation of A e (\)u at the 
boundary of the computational domain we assume u — outside the domain. 
This is legitimate as we choose a large enough computational domain so that the 
solution is essentially zero at the boundary. 

The effective equation (1.5) is solved via a second order centered finite differ- 
ence scheme. For the second derivatives we use the standard stencil (T) 2 w)j := 
(Ax)~ 2 (iUj + i — 2wj + Wj_i) and for the fourth derivatives we use (D i w)j := (Ax) -4 
(wj+2 ~ 4«7j + i + 6wj — 4w; :) _ 1 + Wj_ 2 ) so that the semidiscrete problem in the case 
n — 1 reads 

((I - e 2 ET> 2 )d?u) . = ((AD 2 - e 2 FB 4 )u). . 

We recall that E and F are scalars when n = 1. Discretization in time is performed 
analogously to the case of equation (1.1). 

The above described methods generalize to n > 2 dimensions in a natural way, 
see [6] for equation (1.1) with n = 2. 

In general the parameters a*, a, and /?, which determine the coefficients A, E 
and F in the effective equation, need to be computed numerically. They can be 
computed by numerically differentiating the eigenvalue /io as defined in (2.19). 

4.1 One space dimension 

We choose the material function ay{y) = 1.5 + 1.4cos(y) and the initial data 
f(x) = e~ 0Ax and numerically investigate the quality of the approximation given 
by the effective equation. For the coefficients A = a* and C = a we find 

a* w 0.5385, a w -0.5853, 

so that AD 2 = a*d 2 x w 0.5385 d 2 , ED 2 = -±Cd 2 x w 1.0869^. 

Equation (1.1) was solved with Ax = 27re/30 and At = 0.008 and (1.5) was 
solved with Ax ~ 27r/100 and At = 0.005. In Fig. 1 we plot u e and w £ for 
e = 0.05 at t = 400 = e~ 2 and for e = 0.1 at t = 200 = 2e~ 2 . We see that 
in both cases the main peak and the first few dispersive oscillations are well 
approximated by the effective model. In the latter case, i.e. with t relatively 
large for a given e, a slight disagreement in the wavelength of the tail oscillations 
is visible. Fig. 1 additionally shows oscillations traveling faster than the main 
pulse. These oscillations are physically meaningful as their speed is below the 
maximal allowed propagation speed c := \Y\ J R a Y 1 ^ 2 {y)dy, see [17], marked by 
the vertical dotted line. 

In Fig. 2 we study the convergence of the L 2 (M)— error for the same material 
function and initial data as above. The error is computed at e = 0.2, 0.1 and 0.05 
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Figure 1: One-dimensional equation: the solutions u e and w £ for ay(y) = 1.5+1.4 cos(y) 
and f(x) = e~ 0Ax are compared. Only the right propagating part of the solution is 
plotted. In (a) we have e = 0.05 and in (b) e = 0.1. The insets zoom in on the dispersive 
oscillations to the left of the main peak. 



and t = e~ 2 . The error values are approximately 0.1954,0.0977,0.0494. Clearly, 
the numerical convergence is close to linear, in agreement with Theorem 1.1. 

0.2 

'lu 0.15 

II 

3 0.1 

I 

Is, 

= 0.05 

0.05 0.1 0.15 0.2 

£ 

Figure 2: Convergence of the L 2 -error \\u e — w e \\ L 2 at t = e~ 2 for oy(y) = 1.5+1.4 cos(y), 
f(x) = e~ 0Ax2 , and the three values e = 0.2, e = 0.1, and e = 0.05. We emphasize that 
this is a severe test for convergence: in both steps, e is halved and the time instance is 
quadrupled. 




4.2 Two space dimensions 

Full two-dimensional (n = 2) simulations for small values of e > and time inter- 
vals of order 0(e~ 2 ) are computationally expensive due to the need to discretize 



24 



T. Dohnal, A. Lamacz, B. Schweizer 



each period of size 0(e) x 0(e) in a domain of size 0(e~ 2 ) x 0(e~ 2 ). We therefore 
perform instead a simulation that is designed to mimic the long time behavior of 
a solution originating from localized initial data. After a long time the solution 
develops a large, close to circular, front. Within the strip 

O s := i 6 1 x (— eir, eir) 

we can expect that the front is nearly periodic in the £ 2 — direction. Therefore, we 
perform tests on Q s with periodic boundary conditions in x 2 , and initial data that 
are localized in x\ and constant in x 2 - Our choice is to take f(x) = e~ 0SXl , x e fl s . 
We select a material function that describes a smoothed square structure, namely 

a Y (y)=l + c(y)-c, (4.3) 
1 2 

<y) =g II t 1 + tanh WW + iM)] t 1 ' tanh ( 4 (% - H)] > 

where c := t^t J Y c(y)dy. This choice ensures a relatively large value of the dis- 
persive coefficient a. We find 

a* « 0.5808, a « -0.3078, /3 « 0.0515. 

These values correspond to case 2 in Remark 3.1 so that AD 2 = a* A « 0.5808 A, 
ED 2 = ^A « 0.5300 A, FD 4 = (\a\ + 3/3)8^8^ w 0.4623 Sg^. Due to the 
x 2 — independence of the initial data, the solution of the effective model (1.5) on 
Q s stays constant in £ 2 so that FD A can be dropped and (1.5) becomes 

8 2 w e = 0.5808 8 2 xi w £ + e 2 0.53 d%$w e . 

In the simulations of (1.1) we use Ax\ = Ax 2 = 27re/30 and At = 0.004, and in 
(1.5) we use Ax x = 2tt/100 and At = 0.01. 

In Fig. 3 the main part of the right propagating half of the solution u e is 
plotted for e = 0.1 at t — 100 = e~ 2 . One clearly sees dispersive oscillations 
behind the main pulse. Fig. 4 shows the agreement between w £ and the x 2 — mean 
of u £ at e = 0.1 and t = e~ 2 . 

Conclusions 

We have performed an analysis of wave propagation in multi-dimensional hetero- 
geneous media (periodic with length-scale e > 0). It is well-known that for large 
times, solutions cannot be approximated well by the homogenized wave equation. 
We have provided here a suitable well-posed dispersive wave equation of fourth 
order that describes the original solution u e on time intervals of order 0(e~ 2 ). 
Our analytical results provide an error estimate of order 0(e) between u e and the 
solution w e of the dispersive equation. The coefficients of the effective equation 
are computable from the dispersion relation, which, in turn, is given by eigenval- 
ues of a cell-problem. The qualitative agreement between u e and w e is confirmed 
by one-dimensional numerical tests, that even provide a confirmation of the linear 
convergence of the error in e. In two space dimensions we can observe the validity 
of the dispersive equation in a simplified setting, computing solutions on a long 
strip. 
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Figure 3: Two-dimensional equation: (a) The periodic structure a e {x) given by (4.3) 
over a section of the strip Q s . (b) The main part of the right propagating part of the 
solution u £ at t = 100 for e = 0.1 and f{x) = e -0,61 !. (c) The X2-profile of u e at x\ = x\ 
with x\ being the position of the peak of the pulse. 
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